;;--------------------------------------------------------------------
;; calculate the decaying process of radon and its daughters 
;; exponential decay, reach equilibrium after 40 mins 
;;--------------------------------------------------------------------

   dt   =  1.         ;; time step 1s 
   nd   =  4          ;; number of day simulation  
   nsd  = 24* 3600    ;; second of a day 
   nsh  = 3600        ;; second of a hour 
   nsm  = 60          ;; second of a min
   nsy  = nsd * 365   ;; second of a min
   ndt  = floattointeger(20.*nsd/1.)     ;; number of integration timestep 

   t0   =  5.50 * nsd  ;; efolding time of radon  
   t1   =  4.40 * nsm  ;; efolding time of 218po
   t2   = 38.70 * nsm  ;; efolding time of 214pb
   t3   = 28.40 * nsm  ;; efolding time of 214bi
   t4   =  2.00 * nsd  ;; efolding time of some tracer 

   w0   = 1./t0
   w1   = 1./t1
   w2   = 1./t2
   w3   = 1./t3
   w4   = 1./t4

   v10  = w1/w0 
   v21  = w2/w1 
   v20  = w2/w0 
   v32  = w3/w2 
   v31  = w3/w1 
   v30  = w3/w0 
   v43  = w4/w3 
   v42  = w4/w2 
   v41  = w4/w1 
   v40  = w4/w0 

   c0   = 10.        ;; initial concentration 10 Bq/m3 STP 
   c1   = 0. 
   c2   = 0. 
   c3   = 0. 
   c4   = 0. 

   xa   = new((/ndt/),"float") 
   xb   = new((/ndt/),"float") 
   xc   = new((/ndt/),"float") 
   xd   = new((/ndt/),"float") 

   ;; radioactive decay 

   do it = 0, ndt-1

   c0 = 10. 

   ;; its own decay 
 
   h0 = exp(-dt/t0) * c0
   h1 = exp(-dt/t1) * c1   
   h2 = exp(-dt/t2) * c2 
   h3 = exp(-dt/t3) * c3 
   h4 = exp(-dt/t4) * c4 

   ;; from its parent 

   p0 = 0.               ;; nothing for radon  
   p1 = (c0 - h0)*v10    ;; radon decay to 218po  
   p2 = (c1 - h1)*v21 
   p3 = (c2 - h2)*v32
   p4 = (c3 - h3)*v43 

   ;; total change 
  
   c0 = h0 + p0 
   c1 = h1 + p1 
   c2 = h2 + p2
   c3 = h3 + p3 
   c4 = h4 + p4 

   print("it,c0,c1,c2,c3,c4 :"+it+" "+c0+" "+c1+" "+c2+" "+c3+" "+c4)
 
   end do 

   print("c1/c0,c2/c0,c3/c0 :"+c1/c0+" "+c2/c0+" "+c3/c0)
   print("t1/t0,t2/t0,t3/t0 :"+t1/t0+" "+t2/t0+" "+t3/t0) 
   






